arXiv: 1502.06055v2 [quant-ph] 17Jul2015 


Synchronization of Interacting Quantum Dipoles 

B. Zhu 1 , J. Schachenmayer 1 , M. Xu 1 , F. Herrera 2 , J. G. Restrepo 3 , M. J. Holland 1 ,* and A. M. Rey ^ 
1 JILA, NIST, Department of Physics, University of Colorado, 440 UCB, Boulder, CO 80309, USA 
2 Department of Physics, Universidad de Santiago de Chile, 

USACH, Casilla 307 Correo 2 Santiago, Chile and 
3 Department of Applied Mathematics, University of Colorado, Boulder, Colorado, 80309, USA 

(Dated: July 20, 2015) 

Macroscopic ensembles of radiating dipoles are ubiquitous in the physical and natural sciences. 

In the classical limit the dipoles can be described as damped-driven oscillators, which are able to 
spontaneously synchronize and collectively lock their phases in the presence of nonlinear coupling. 

Here we investigate the corresponding phenomenon with arrays of quantized two-level systems cou¬ 
pled via long-range and anisotropic dipolar interactions. Our calculations demonstrate that by 
incoherently driving dense packed arrays of strongly interacting dipoles, the dipoles can overcome 
the decoherence induced by quantum fluctuations and inhomogeneous coupling and reach a synchro¬ 
nized steady-state characterized by a macroscopic phase coherence. This steady-state bears much 
similarity to that observed in classical systems, and yet also exhibits genuine quantum properties 
such as quantum correlations and quantum phase diffusion (reminiscent of lasing). Our predictions 
could be relevant for the development of better atomic clocks and a variety of noise tolerant quantum 
devices. 

PACS numbers: 05.45.Xt, 42.50.Lc, 37.10.Jk, 64.60.Ht 


I. INTRODUCTION 

Arrays of synchronized oscillators [1] are ubiquitous in 
biological [2, 3], physical [4] and engineering [5] systems 
and are a resource for technological advances [6]. Al¬ 
though there has been significant progress in the study of 
synchronization in classical systems [7], the understand¬ 
ing of the same phenomena in the quantum realm re¬ 
mains limited. A major obstacle so far is the general 
problem of the exponential scaling of the Hilbert space 
with system size which makes calculations dealing with 
quantum arrays very challenging. In fact, current investi¬ 
gations have been limited to the exact treatment of arrays 
of a small number of coupled quantum oscillators [8-18], 
and large ensembles at the mean field level or by includ¬ 
ing quantum corrections perturbatively [19-21]. Highly 
symmetric situations with collective coupling mediated, 
for example, by a cavity mode [22-24], have also been 
studied. 

Ensembles of radiating dipoles are a natural platform 
to study quantum synchronization, where coherence can 
be generated from an incoherent source. One might re¬ 
gard laser systems, where radiation is amplified by the 
stimulated emission of photons, as a prototypical exam¬ 
ple. However, lasing is fundamentally a distinct phe¬ 
nomenon from quantum synchronization. This can be 
seen from the fact that lasing is possible even in the ab¬ 
sence of coupling between the atomic dipoles, as is clear 
in the single atom laser [25], or in atomic beam lasers 
where only one atom is present in the cavity at any given 
time. A more relevant situation is the quantum syn¬ 
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chronization that takes place in the context of superradi¬ 
ance [26]. It has recently been understood that, in con¬ 
trast to lasers, steady-state superradiance can produce 
spectrally pure light [26-28] without stimulated emission. 
So far this has been demonstrated using a cavity mode as 
a communication channel that spatially selects an optical 
mode and enhances the coupling (through the cavity fi¬ 
nesse). A more generic and relevant scenario, with great 
potential and applicability, is the emergence of sponta¬ 
neous macroscopic quantum synchronization in radiating 
dipole arrays without a cavity but naturally coupled by 
the intrinsic anisotropic and long-range dipolar interac¬ 
tions. This is the situation considered in this paper. 

Here we demonstrate that in the presence of an in¬ 
coherent repumping source, dipole induced cooperative 
emission can dominate over spatial inhomogeneities and 
quantum fluctuations and lead to a resilient steady-state 
that exhibits macroscopic quantum phase coherence and 
intrinsic quantum correlations. An iconic example of a 
macroscopic coherent state is a Bose-Einstein conden¬ 
sate, achieved in ultra-cold gases at thermal equilibrium. 
In our case, however, the macroscopic order is reached in 
the steady-state of an interacting and driven-dissipative 
system. Moreover, the cooperative behavior can be de¬ 
tected by measuring the spectral purity of the emitted 
radiation. We note that in clear distinction to previous 
studies [8, 15-17, 19, 21], our proposal does not rely 
on an external coherent source or externally generated 
nonlinearities to seed the collective phase. In our model 
synchronization emerges as a spontaneously broken sym¬ 
metry driven by incoherent processes in naturally coupled 
dipole arrays. As we show, and somewhat counterintu¬ 
itively, an incoherent drive is sufficient to generate phase 
coherence in these systems. 

Specifically, the systems we consider are dense arrays 
of frozen quantum dipoles modeled as quantized two-level 
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(a) Cold atoms in a crystal of light (b) Inelastic interaction Elastic interaction 



FIG. 1. Arrays of quantum dipoles spontaneously emit and absorb photons at rate T. The photons mediate dipolar interactions 
between dipoles separated by a distance r a b with both dissipative, /( r a &), and elastic, g( r a &), components. A repumping source 
at a rate W provides energy to maintain the oscillations and can be implemented using additional internal states that are 
not shown, (a) A possible implementation using cold atoms in an optical lattice, (b) The dipolar couplings g(r a b) and f(r a b) 
exhibit a complex angular distribution as a function of 6, the angle between the dipole orientation (determined by an external 
electromagnetic field) and |r a &|. The maximum value of / and g for fixed |r a b| is denoted as /max and g max . The cone illustrates 
the magic angle, 6m — arccos(l/v / 3)- 


systems. By dense arrays of frozen dipoles we mean ar¬ 
rays separated by a distance much closer than the wave¬ 
length of the emitted photons and with motional de¬ 
grees of freedom evolving at a much slower rate than 
their internal dynamics (Fig. 1). These conditions can be 
readily satisfied in a variety of quantum systems found 
in atomic, molecular and optical physics (e.g., Rydberg 
gases [29-31], alkali vapors [32], alkaline-earth atoms [33], 
and polar molecules [34]), chemistry (e.g., J-aggregates of 
dye molecules [35-37]), and biology (e.g., light-harvesting 
complexes [38, 39]). In cold vapors, one possible way to 
freeze the motion and tightly trap the particles is via 
an optical lattice potential (Fig. 1). In this case a sub- 
optical-wavelength transition must be used in order to 
reach the tight-packing regime [33, 34]. 

To fully understand synchronization in the complex 
dipolar system, we analyze each of the ingredients that 
compete and affect synchronization in a step-by-step pro¬ 
cedure: the interplay between repumping and collective 
emission, inhomogeneity in the coupling constants, quan¬ 
tum correlations, and the competition between elastic 
and inelastic interactions. The paper is organized as 
follows: In Sec. II we introduce the system in consid¬ 
eration and the master equation we use to describe the 
dynamics. In Sec. Ill we first provide a simple mean- 
field description and discuss connections to the classical 
Kuramoto model—the iconic model used to describe syn¬ 
chronization in non-linear coupled oscillators. In Sec. IV 
we discuss the phase diagram for the quantum system 
assuming collective (all-to-all) coupling and compare it 
with the mean-field solution. For this exactly solvable 
case we are able to explicitly quantify the entanglement 


and correlations present in the steady state. In Sec. V 
we study how inhomogeneity in the inelastic couplings af¬ 
fects synchronization and focus on the case of power-law 
decaying interactions. In Sec. VI we study the emergence 
of quantum synchronization in radiating dipoles taking 
the full long-range and anisotropic dipolar interactions 
into account. In Sec. VII we discuss experimental imple¬ 
mentations of our model, and in Sec. VIII we provide a 
conclusion and an outlook. 


II. DIPOLE-DIPOLE INTERACTION AND 
MASTER EQUATION 

In this work we consider arrays of quantum dipoles 
with two accessible levels, which we denote as |^) and |/). 
The interactions between two dipoles a and b are de¬ 
scribed by the functions g(r a b) and /(r a &), which de¬ 
pend on the dipoles’ separation, |r a &|, and the angle 0 
between the mean dipole moment and the vector joining 
the dipoles [See Fig. 1(a)] [40]: 

/ \ 3r f . 2 COsCab . to 2 q -.NfCOsCab , sinCafc, \ 
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Here, ( a b = 27r|r a ^|/A, where A is the characteristic wave¬ 
length of the dipole-transition, and F = /(0) is the spon¬ 
taneous photon emission rate from a single dipole. The 
function g(r a b) describes the elastic dipole-dipole interac¬ 
tions, while /(r a fr) gives rise to inelastic collective photon 
emission. These terms are similar to those that determine 
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the radiation of classical electric dipoles, and the depen¬ 
dence on |r a &| reflects the propagation of photons from 
one atom to another. The terms oc 1 /( ab account for re¬ 
tardation effects in the far-field regime and those oc 1 / 
account for instantaneous propagation in the near-field. 
When ( ab < 1, the elastic g interactions with a strong 
angular variation are dominant except close to the magic 
angle 0 m = arccos(l/v / 3), at which they are greatly sup¬ 
pressed. In contrast, /( r ab ) is almost isotropic in the 
near-field regime [see Fig. 1(b)]. 

The spatially uniform behavior of /( r ab ) at short dis¬ 
tance is what gives rise to cooperative effects and super¬ 
radiant emission [26]. Under generic conditions, however, 
superradiance is a transient effect that substantially lim¬ 
its the lifetime of dipole excitations. To compensate for 
the fast decay here we add an incoherent repumping driv¬ 
ing term at a rate W. This term is needed to generate 
a synchronized steady state where long-lasting coherence 
persists. An incoherent repumping drive is commonly 
used in laser systems to maintain population inversion. 
It can be implemented by coherently driving, at a rate 
U ex , the 1 1) state to an excited level that spontaneously 
decays, at a rate 7 U ex , to the state |t)- Due to 
the fast depletion of the excited state, it can be adia- 
batically eliminated and thus the net process is just an 
incoherent transfer of population from ||) to |t) at a rate 

w = n 2 ex /j [4i]. 

The evolution of N dipoles is modeled by a quantum 
master equation for the reduced density matrix p of the 
dipoles [26]: 


f t ~~ l -[H 0 ,p] + C f [p\ + C w [p], 


h N N 

[ S ^a + J2 9(rab)^t^b 

b=l,bj^a 


a= 1 


(1) 

( 2 ) 


C Ap\ = oE /( r ab)(2<3- b pert -°t°b P- P & t°b )> (3) 


a,b 


£w[p\ = T ^( 2 <x+ per a -a a &tp- per a d +). 


( 4 ) 


The Hamiltonian Hq generates the coherent evolution of 
the dipole array where da ’ ’ ’ are the Pauli spin op¬ 
erators for dipole a, S a denotes its bare oscillation fre¬ 
quency and h is the reduced Planck constant. The Lind- 
blad operator functionals, Cf,w, describe the inelastic 
photon emission and incoherent repumping processes, re¬ 
spectively. 


III. MEAN-FIELD TREATMENT AND 
CONNECTION TO THE KURAMOTO MODEL 


To obtain a qualitative picture of how synchroniza¬ 
tion can happen among the dipoles, we first perform 
a mean-field treatment and show the close connection 
between our quantum model and the prototype models 


for classical synchronization. The mean-field approach 
assumes uncorrelated dipoles, i.e., p = ® a /U, where 
each p a = ^crcr / =t iPa ,a l cr )( (j/ | a 2 x 2 matrix in 
the pseudospin 1/2 basis {|t) ? |l)}- The components of 
the single-dipole density matrix, p a , can be visualized as 
a Bloch vector {S^{t) cos0 a (t), S^{t) sin0 a (t), S%(i)} = 
(4/2){/7 + pit, -i(pi l - pi t ),pi t - pi 4 '} (Fig. 2). The 
mean-field solution yields a system of coupled nonlinear 
differential equations for p° ,a . For each a = 1,2,..., AT 
the parameters evolve as 

rjQ z (f\ _ r 

a =e'-St(t) yt S b(t) f(r a b)cos[5<t) b a\-g(rab)sm[5<t) b a 


dt 


-r 


b^a 


2 +S Z a(t) 


+ w 


2 -^) 


( 5 ) 


dSt(t) 


dt 


s Z( t )’V) S b( t ) g( r ab)sm[5(/) ba ]- f(r ab )cos[8(j) ba ] 


b^a 


r + w 




(6) 


dd>a 


N 




b=l,b^a a 


g(r ab )cos[S(/) ba } -\-f(r ab )sm[S^ ba 


( 7 ) 


where S(p ba (t) = <fi b (t) — 4> a (t). The term proportional 
to f(r ab ) in Eq. (7) that contains the sine function can 
be identified with a similar term in the Kuramoto model 
(KM) [42]: 

7 / N 

= S a + K ^2 sin[<ty ba ], ( 8 ) 

b= 1 

where IT, the coupling strength per oscillator, must be 
large enough and positive for synchronization to oc¬ 
cur. The term proportional to g(r ab ) that contains 
the cosine function appears in the Sakaguchi-Kuramoto 
model [43]— a more general but similar synchronization 
model to the KM. Compared to the basic KM, the situa¬ 
tion here is more complex. This is due to the fact that in 
Eq. (7) the coupling constants are nonuniform and effec¬ 
tively time-dependent, since S^(t) and S*(t) are dynamic 
variables. 

To investigate whether the mean-field model admits 
spontaneous synchronization we consider first the simpli¬ 
fied case where S a = 0 for all dipoles, impose g(r ab ) = 0 
for all pairs, and assume a constant collective decay rate 
Nf(r a b) = / e ff- We define a global order parameter Z as 
Ze l ® = S'^ L e^ a and look for a solution in which Z 

is time-independent and synchronized oscillators possess 
a collective frequency cJ, and thus a macroscopic phase 
4> = uJt. These conditions lead to two equations for the 
order parameter Z and the collective frequency c 0 (see 
Appendix A): 


uj = 0 , 

y/f eB (w-r)-(w + r)* 

\/2/eff 


( 9 ) 

(10) 
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FIG. 2. (a) Mean-field phase diagram calculated from the order parameter 0 < Z < 1/^/8- The insets show snapshots of 
the tips of the Bloch vectors (red points) for dipoles prepared with random initial phases and then evolved to steady-state in 
both regimes, (b) Quantum phase diagram calculated from 0 < Zq < l/\/8. (c) The time evolution of the conditional QFI 
exhibits entanglement (dashed line: single trajectory, solid line: mean value of a few trajectories). Panel (d) shows the steady 
state QFI vs W/T after averaging over many trajectories. The solid line corresponds to the conditional case, and indicates 
entanglement over the repumping range where synchronization exists. Upon computing the ensemble average one recovers the 
reduced density matrix which leads to a calculated QFI below the entanglement witness threshold (dashed line), (c) and (d) 
are shown for / e ff = 15T. For all panels, S a — gi?ab) — 0 and /(r a &) = f e s/N. 


The solution is shown in Fig. 2(a). The insets show 
the phase distribution in the steady-state for an array 
of oscillators initially prepared with a random distribu¬ 
tion of phases for two different values of the repump rate 
W. For a slow repumping rate (bottom inset), the sys¬ 
tem remains unsynchronized. As the repumping rate is 
increased beyond a threshold value, the system enters a 
synchronized state, as can be seen by the appearance of 
phase locking and the resulting narrow phase spread (top 
inset). This can be explained by the fact that one neces¬ 
sary condition for synchronization in the KM is K >0, 
which translates to the requirement S z > 0 on average 
in our model and thus the need to have sufficiently large 
repump rate. In the limit f e ff T (e.g., for large N), 
maximum synchronization is achieved at W op t = / e ff/ 2 , 
where the order parameter Z reaches a maximum value 
Zmax ~ \j 1/8. For this optimal condition for synchro¬ 


nization the quantum dipoles are ordered with the same 
phase and radiate with atomic inversion S z a ^ 1/4. Note 
that the maximum order parameter is smaller than 1/2 
even when fully synchronized because of this required fi¬ 
nite value of the atomic inversion. One intriguing aspect 
is that repumping, which is the process that builds up 
synchronization, is itself an incoherent process. It is cru¬ 
cial that repumping does not preserve the norm of the 
collective Bloch vector, allowing it to extend or contract. 
For large W > W Qpt , Z decreases again reflecting a sup¬ 
pression of synchronization. In this limit the dipoles are 
repumped so fast that they are all driven to the |t) state 
(S z —» 1/2 and —)> 0 ) and phase coherence between 

them cannot build up. 

The cases of a heterogeneous distribution of £ a ’s or fi¬ 
nite g( r a fr) 7 ^ 0 can be treated at the mean-field level in 
a simple way. The results, summarized in the Appendix 
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A, are qualitatively similar. In this case we define the 
effective couplings as Nf(r ab ) = / eff and Ng(r ab ) = g eff . 
In general, the inclusion of a finite spread A in S a de¬ 
creases the value of the order parameter Z. For in¬ 
stance, if S a is sampled from a Lorentzian distribution 
p(S a )=A/[7r(A 2 + S 2 a )}, 


y/feaP -Q 2 + 2A 2 - 2A^A2 + / eff P 


( 11 ) 


where Q = Y + W and P = W — Y. Optimal syn¬ 
chronization is obtained at a smaller repumping rate, 
W op t ~ /eff/ 2 - A/y/2. We note that for given f eS 
and IF, synchronization is destroyed (that is, Z = 0) 
at A c = ( Q 2 - / e ffP)/( 2 Q). 

When the dipoles have identical detunnings, S a = 0 , 
the elastic couplings simply induce a global frequency 
shift cJ = g e ffQ/(2/ e ff) that can be eliminated by moving 
to a rotating frame. 


IV. QUANTUM SYNCHRONIZATION FOR 
THE COLLECTIVE SYSTEM 

In the simplified case where S a = 0 for all dipoles, 
g(r ab ) = 0 fc> r all pairs, and a constant collective decay 
rate Nf(r ab ) = / e fr, it is possible to exactly solve Eq. ( 1 ), 
i.e., the full quantum dynamics, even for many particles, 
allowing us to benchmark the validity of the mean-field 
solution. This is due to the invariance of the master 
equation under individual dipole permutations that re¬ 
duces the scaling of the Liouville space from exponential, 
4^, to polynomial, of order N 3 [44]. 


A. Phase Diagram 

Quantum fluctuations can lead to phase diffusion and 
to decay of single particle coherences in the steady-state 
(it is possible for (<r+) —» 0 even in a synchronized state), 
so Z cannot be used as a measure of synchronization in 
a beyond mean-field treatment. However, phase locking 
in quantum mechanics can be quantified by the degree 
of spin-spin correlations Zq, defined by Zq = (crt(J b ), 
where the bar indicates an average over all pairs of dif¬ 
ferent dipoles a and b. For an unsynchronized state 
Zq is 0 and for a completely synchronized state Zq is 
Zg ax = 1/y/S [27, 28]. The corresponding phase dia¬ 
gram, shown in Fig. 2(b), closely resembles the mean- 
field one. 

To demonstrate that Zq can be used to quantify the 
emergence of quantum synchronization, regardless of the 
inherent non-equilibrium and dissipative character of our 
system, we have also computed pairwise two-time correla¬ 
tion functions (see Appendix B). The decay rate of these 
correlations encodes information about the spectral co¬ 
herence of the emitted radiation. The range of W/Y val¬ 
ues where the emitted light is maximally coherent agrees 


with the regime where the system is optimally synchro¬ 
nized according to Zq. Moreover, we have also confirmed 
the moderate importance of higher order correlations in 
the synchronized steady-state by comparing the exact so¬ 
lution with a cumulant expansion calculation (which in¬ 
cludes lowest order corrections to the mean-field result). 
We find the cumulant expansion agrees well with the 
exact solution (see Appendix B). The only limit where 
there are important deviations is at very weak pump¬ 
ing IF T where the system favors subradiant emission 
arising from strong atom-atom correlations [indicated by 
the purple region in Fig. 2(b)] [44]. 


B. Quantum correlations and entanglement 

The robust macroscopic quantum coherence exhibited 
by the synchronized state leads to the natural question 
of whether or not entanglement can be present in the 
steady-state even in this dissipative environment. Most 
previous studies that attempted to address this question 
have been limited to small systems [13, 18, 21, 45] and 
focused on the entanglement between a pair of synchro¬ 
nized oscillators. Here, to determine the non-separability 
of the many-body steady-state, we compute the average 
of the quantum Fisher information (QFI) and use it as an 
entanglement witness [40, 46, 47]. Any N particle state 
with ( N 2 + 2N)/3 > Fq(p) > 27V/3 is entangled (non- 
separable) and a quantum resource for phase estimation 
(see Appendix C for details). 

Due to dissipation, the density matrix of the system 
is reduced to a mixed state, as obtained from Eq. (1), 
which describes the dynamics of the system after a sta¬ 
tistical average over many experimental trials. However, 
the evolution of the system for an individual experi¬ 
mental realization can be quite different. We consider 
a Gedanken experiment in which one monitors the sys¬ 
tem evolution and keeps a measurement record of the 
emitted photons. The evolution of the system is then 
conditioned on the measurement record [48, 49]. This 
type of conditional evolution has been widely studied in 
quantum optics and utilized together with quantum feed¬ 
back control in examples such as the optimal generation 
of spin squeezed states (See Refs. [50, 51]). It should 
be emphasized that the conditional evolution based on 
the measurement record gives a quantum trajectory that 
should not be regarded simply as a numerical tool to al¬ 
low the efficient assembly of ensemble averages. Each 
quantum trajectory is a potentially realizable physical 
outcome (even if hard to perform in practice) as allowed 
by the quantum dynamics of the open quantum system 
under consideration. 

We calculate Fq(p c ) (with the c in p c meaning condi¬ 
tional) for each conditional trajectory and in Fig. 2(d) 
we show its average over a sufficiently large set of trajec¬ 
tories at steady state (see Appendix C). For this condi¬ 
tional evolution we observe entanglement in a parameter 
regime that correlates with Zq > 0 [see Fig. 2(c) and (d)]. 
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FIG. 3. Quantum correlations and total correlations. The total and quantum correlations in the steady state are quantified 
by the mutual information 0 < X < 2 and quantum discord 0 < T> < 2. (a) In the synchronized phase, there are nonzero 
quantum correlations and classical correlations (X — V) , and both show a dependence on W that qualitatively agree with Zq . 
(b) Even in the thermodynamic limit, quantum correlations remain a significant fraction of the total correlations. For both 
panels, S a = g(r ab ) — 0 and /( r ab ) = f e s/N. 


On the other hand, if we discard the information present 
in the measurement record, by using the ensemble aver¬ 
aged p obtained from directly solving Eq. (1), and then 
computing Fq(p ), the QFI falls below the entanglement 
witness threshold [see Fig. 2(d)]. 

To differentiate quantum effects from classical ones, we 
further calculate the quantum discord V [40], which can 
be considered as a measure of quantum correlation more 
general than entanglement and more robust in a dissipa¬ 
tive environment [52-54]. Separable states with nonzero 
V are intrinsically nonclassical, since local measurements 
performed on a subsystem inevitably disturb the whole 
system [52, 55]. We measure classical correlations of the 
steady state by the difference between the mutual infor¬ 
mation X [40] and V. We find the mixed steady-state con¬ 
tains nonzero quantum correlations in the synchronized 
regime Fig. 3(a). Moreover we observed that although 
both V and Zq exhibit a similar dependence with pump¬ 
ing rate W, they do not exactly peak at the same value 

[15]- 

Although the existence of a nonzero V has been re¬ 
ported to exist in several quantum synchronization stud¬ 
ies [18, 55], we want to emphasize that those have been 
always limited to small systems. To our knowledge our 
calculations are the first to consider V in macroscopic 
samples. In Fig. 3(b) we show the dependence of V with 
system size. Our calculation shows that quantum cor¬ 
relations remain a significant fraction of X even in the 
thermodynamic limit. 


V. SYNCHRONIZATION WITH 
FINITE-RANGE INTERACTIONS 

Up to this point we have only considered all-to-all in¬ 
teractions; now we consider the effect of finite range in¬ 
teractions on synchronization. In the dipole array both 
f(r a b) and g(r a fr) are nontrivial functions of r a ^ and 


contain terms decaying as a power-law with distance, 
oc l/|r a fr| a with a = 1,2,3. Instead of dealing with all 
these terms together, to gain insight on how spatial in¬ 
homogeneities affect quantum synchronization, we first 
study a simpler case assuming a power-law cooperative 
decay /(r a &) oc |r a ^| _a , with the exponent a as a variable 
parameter and set both g( r a ^) = 0 and 5 a = 0. 

In the classical regime [56-59] , analytical calculations 
and numerical simulations considering arrays of oscil¬ 
lators interacting via power law interactions on a one¬ 
dimensional lattice had identified a c = 3/2 as the critical 
value of the power law exponent below which long-range 
phase order is possible [58]. For a < <a c , a transition to 
a state in which a finite fraction of the oscillators is en¬ 
trained takes place for a sufficiently strong but finite cou¬ 
pling strength in the large system limit. Generalizations 
of these results to oscillators in D dimensions [58] have 
also identified three different regimes for synchronization: 
perfect phase ordering for a < D , entrainment with long- 
range phase order for a < 3D /2 and a crossover to expo¬ 
nential decay of correlations at a = (3D+l)/2. Reference 
[59] has also suggested that in the regime a > D global 
synchronization is absent but local synchronization per¬ 
sists for arbitrary weak coupling with a slowly decaying 
order parameter. 

To quantify the effect of finite-range interactions on 
synchronization in the quantum regime we compute spin- 
spin correlations within linear clusters that contain d 
dipoles, (Zq) 2 = (frad'b)d, using a cumulant expansion 
method as described in Appendix B. Here the bar fol¬ 
lowed by a subscript d indicates an average over the pairs 
of different dipoles a and b contained in a linear cluster 
of size d. The linear clusters start at the central spin as 
shown in Fig. 4. We have confirmed that the cumulant 
expansion method reproduces well the correlation func¬ 
tions by performing direct comparisons with the exact 
solution (see Appendix B). Fig. 4(a) shows the behavior 
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FIG. 4. (a) Spin-spin correlations, Zq, in linear clusters containing d dipoles at optimal repumping W for power-law couplings 

f(r a b) — j ( : Ff) lattice spacing a. We set S a = g{?ab) — 0 and consider N — 900 dipoles arranged in both linear (D = 1) 

and square lattice (D — 2) geometries. For a < D global synchronization is observed and the order parameter is independent 
on cluster size d. For D < a the order parameter starts to clearly decay with increasing d. The magenta line (Zq = 0.14) 
provides an indicative scale of the boundary between global and local synchronization. The white contour lines provide an 
indication of the decrease of the synchronized domains with increasing a. (b) Pair wise two-time correlation functions in the 
steady-state are parametrized by Z a ^(r) — Acosiyr) exp(— jr) where a is chosen as the central dipole of a linear chain of 
N = 200 dipoles. The dipoles are assigned random detunings S a distributed uniformly in [—T/2,r/2]. The dark blue, red, and 
light blue symbols correspond to a = 0, 0.65 and 2 respectively. The histogram of frequencies v exhibits similar synchronization 
regimes than those seen in (a). 


of Zq as a function of cluster size d and power-law decay 
exponent a in arrays of dimension D = 1 and 2 . Clear 
global synchronization with an order parameter Zq in¬ 
dependent of d is observed for a < D. For 0 < a < D/2, 
the local order parameter Zq is almost independent of 
a and d and the system behaves almost like the all-to- 
all system. For D/2 < a < D synchronization remains 
global and almost independent of d, but the order pa¬ 
rameter slowly decreases with a. For a > D, synchro¬ 
nization becomes local and correlations quickly decrease 
with cluster size. The magenta contour provides an in¬ 
dicative scale of the boundary between global and local 
synchronization. The white contour lines also provide in¬ 
formation about the decrease of the order parameter with 
increasing a and d. We observe that, as in the classical 
case, a ~ D roughly marks the transition between global 
and local synchronization, although a more quantitative 
comparison would require far larger systems. 

An alternative way to characterize domain formation 
and the fact that it can persist even when there is a 
variation in the local detunings, S a ^ 0 , is to exam¬ 
ine pairwise two-time correlation functions, Z a ^(r) = 
lim t ^ 00 ((a+(t + t) + (t + T))(d~ (t) + a^(t))), which 
can be related to the emission spectrum of the pair of 
atoms [48]. The oscillations in Z a ^(r) encode infor¬ 
mation about the relative precession rate between dif¬ 
ferent dipoles. By parameterizing Z a ^(r) as Z a ^(r) = 
Acos(z/ a 5 r) exp(—yr) we can extract the relative preces¬ 
sion frequency v^ between dipoles a and 5, where en¬ 
trainment of dipoles a and b corresponds to v a \) = 0 . 
To explore the entrainment of dipole pairs in our sys¬ 
tem, we assign random detunings distributed uniformly 


in [— r/2, r/2] to a linear chain of N = 200 dipoles and 
calculate for b = 1 , 2 ,..., 100 with a = 101 corre¬ 
sponding to the central dipole in the chain. Synchro¬ 
nization regimes similar to those shown in Fig. 4(a) are 
observed for this Dm 1 system, which we illustrate in 
Fig. 4(b) by plotting a histogram (top panel) and the dis¬ 
tribution of frequencies v (bottom panel) for three values 
of a. For global coupling, a = 0 (dark blue symbols), all 
the dipoles become entrained with each other [y = 0 ), 
indicating complete synchronization; for a = 0.65 (red 
symbols) dipoles split into entrained (y = 0 ) and drifting 
(y 7 ^ 0) groups. While not all dipoles are entrained, the 
entrained dipoles are distributed along the whole array, 
and thus synchronization is still global; and for a = 2 
(light blue symbols) the majority of dipoles are not en¬ 
trained. These observations of relative precession fre¬ 
quencies between pairs of oscillators are consistent with 
the regimes obtained from the order parameter plotted 
in Fig. 4(a). 


VI. SYNCHRONIZATION OF DIPOLES WITH 
ELASTIC INTERACTIONS 

We now treat the full problem of radiating quan¬ 
tum dipoles incorporating elastic interactions g(r a fr) and 
the intricate competition of spatially-dependent and 
anisotropic couplings [both g(r a fr) and /(r a &) have terms 
with power law dependence a = 1,2,3 on distance] 
(Fig. 1 ). We solve the full master equation without any 
approximation [48] for systems of up to twenty dipoles 
in a chain using the actual spatial dependence of both 
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FIG. 5. (a) Synchronization in dipole arrays is demonstrated 
for N = 12 dipoles on a line when subjected to incoher¬ 
ent pumping (optimal rate). In this geometry, regardless 
of the strong angular variation of g with the lattice spac¬ 
ing a (see contours) the order parameter, Zq , (normalized by 
Zg ax = I/a/ 8) exhibits a weak dependence on 6 and a and 
reaches a maximum at 6 — 6m- (b) The order parameter is 
computed for N = 16 dipoles on a line with 6 — Om and 
/eff = J2a,b^af( r ^)/(N - 1) (symbols), and for a system 
with constant f(r a b) — fe ff/iV and g(r a b ) = 0 (dashed line)s. 
Similar dependence on W is found for these two different sys¬ 
tems. Here the order parameter for dipoles is always smaller 
in the presence of elastic interactions. 


/(r a b) and g(r a &), and set S a = 0. We observe a robust 
synchronized state t hat exi sts in a wide parameter space. 
As long as / e ff = Nf(r a b) is large enough, we find that 
synchronization takes place and is only weakly affected 
by substantial differences in g( r a &), e.g., variations in the 
dipole array that modify g(r a b) by two orders of magni¬ 
tude only decrease the order parameter by a factor of two 
(Fig. 5) in the steady state. This is in striking contrast to 
the situation in a system without dissipation, where the 
elastic interaction is known to generate entanglement be¬ 
tween spins and to cause a decay of the order parameter 
during time evolution [60]. 

For the orientation = arccos(l/A/3), the order pa¬ 
rameter reaches a significant fraction of ^g ax , indicat¬ 
ing the emergence of macroscopic spontaneous synchro¬ 
nization of the radiating quantum dipole array (Fig. 5). 
To further emphasize the relevant role played by the in¬ 
elastic term, in Fig. 5(b) we compare a solution of the 
master equation [Eq. (1)] for two cases: a system of cou¬ 
pled dipoles arranged in a ID chain and oriented at the 
magic angle (symbols) and an array of identically coupled 
dipoles with the same / e ff but experiencing only inelastic 


interactions \g(r a b) = 0, dashed lines]. The calculated or¬ 
der parameters agree well for the two different cases. The 
similar behavior demonstrates that in spite of the com¬ 
plex geometry of the dipolar interactions, the capability 
of the dipole system to synchronize can be characterized 
to great extent by the quantity f e ff. 


VII. EXPERIMENTAL IMPLEMENTATION 

Our calculations above demonstrate the potential for 
synchronization in a dense array of dipoles. The flexible 
and precise control exhibited by ultracold atomic systems 
make them ideal platforms to experimentally investigate 
the synchronization phenomenon predicted here. Atomic 
systems operate with a large number of quantum oscil¬ 
lators and also allow for the tunability of the interaction 
parameters over a broad range. 

One possible set-up to observe synchronization con¬ 
sists of arrays of ultracold 87 Sr atoms prepared in two 
electronic internal states that form the two-level system. 
The |/) could then correspond to the long-lived 5s5p 3 Pq 
state, with an intercombination line narrower than 10 -3 
s -1 . This is the state used to operate the most pre¬ 
cise atomic clocks [61]. The |t) could correspond to the 
5s4d s Di state with a natural linewith T = 290 x 10 3 
s -1 . Both states can be trapped in an optical lattice 
at the magic wavelength a = 0.2 pm [33], that gener¬ 
ates the same trapping potential for both states mini¬ 
mizing Stark shifts and inhomogeneities in the coupling 
constants. The dipole-dipole interactions are mediated 
by photons at the wavelength A = 2.6 pm and thus, as 
shown in Fig. 5, the ratio a/X < 0.08 falls in the param¬ 
eter regime where dipoles can be synchronized. 

By changing the angle between the laser beams used 
to form the lattice potential, the lattice spacing can be 
varied allowing tunability of the interaction strength be¬ 
tween dipoles. The incoherent pumping can be realized 
by coherently transferring the 5s5p 3 Po population to 
one or several appropriate intermediate states that decay 
rapidly to the 5s4d 3 D i state [62]. An example would be 
the 4d5p 3 Pi state [63]. 

The polarization of the dipoles can be oriented in an 
arbitrary direction by an electromagnetic field. Although 
all the dipoles cannot be oriented at the magic angle in a 
3D geometry, one may still suppress the elastic interac¬ 
tions by dynamical decoupling techniques adopted from 
NMR [64]. Those have been already demonstrated in 
ultracold polar molecule systems [34]. Another possi¬ 
bility is to use a spatial configuration of external fields 
that induces an ‘averaging out’ of the dominant elastic 
interactions [65]. Moreover, by slightly departing from 
the magic-wavelength condition, the dipoles can be sub¬ 
jected to onsite inhomogeneities that generate different 
detunings S a . 

The phase synchronization can be probed by measur¬ 
ing Zq, which experimentally can be directly obtained 
from the fluorescence intensity. As suggested in Sec. V, 
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phase locking can also be extracted from two-point cor¬ 
relations which can be determined by analyzing the flu¬ 
orescence spectrum [62]. 

An intriguing but also more speculative and less con¬ 
trollable realization of our quantum dipole model is the 
case of fluorescent organic molecules. A possible two- 
level configuration in those systems consists of a vibra¬ 
tional level of the ground electronic potential chosen as 
1 1) and the lowest vibronic level of the first excited po¬ 
tential chosen as |t). Incoherent pumping can be real¬ 
ized by driving an optical transition to a higher excited 
vibronic level \<fi) in the first excited potential, which de¬ 
cays on picosecond timescales to the state |t) via non- 
radiative transitions [ 66 ]. Typical values of the fluo¬ 
rescence wavelength Af and lifetime Tf for organic chro- 
mophores under a variety of environmental conditions 
put these systems in a regime of near optimal synchro¬ 
nization. For instance, pseudoisocyanine chloride (PIC) 
and merocyanine derivatives commonly used in organic 
light-emitting diodes (LED) [67-70] typically form low¬ 
dimensional molecular aggregates in liquid solution with 
a ^ 0.5 — 2.0 A, and ratios a/Af on the order of 10 -3 . 
The typical fluoresence decay rate for these organic chro- 
mophores is T ^0.1 — 1 GHz [ 66 ]. In order to achieve 
W/Y = 1 and enter the synchronized phase, the required 
pumping laser intensity is Iw ~ 1 — 10 kW/cm 2 , which 
is lower than the theoretical lasing threshold intensities 
It h - 0.1 - 1 MW/cm 2 of dye lasers [ 66 ]. Therefore, it 
should be feasible to achieve steady state synchronization 
of organic dipoles via incoherent optical driving. 


VIII. CONCLUSION 

We have demonstrated that a system of radiating 
quantum dipoles can be synchronized in the presence 
of repumping. Our analytic mean-field approach pro¬ 
vides a direct analogy between synchronization of quan¬ 
tum dipoles and synchronization of classical phase os¬ 
cillators. Using exact solutions of the master equation 
and a cumulant expansion approach, we determined the 
necessary conditions for synchronization, and the entan¬ 
glement properties in the steady state of macroscopic en¬ 
sembles under different measurement protocols. We also 
analyzed the effect of finite-range interactions in large 
arrays. To our knowledge those have been previously ex¬ 
plored only in the classical regime. For treating the gen¬ 
eral case of dense packed dipoles, we numerically solved 
the master equation exactly for up to twenty dipoles, and 
studied the effect of anisotropic elastic interactions. 

Our results show that the intrinsic macroscopic co¬ 
herence of the superradiant steady state is inherently 
resilient to single particle decoherence, spatial inhomo¬ 
geneities, and noisy environmental effects. This observa¬ 
tion could have relevant application to the development 
of low-threshold organic lasers, highly efficient solar cells, 
materials with enhanced chemical reactivity, as well as 
ultra-precise quantum devices, where these effects are 


anticipated to play an important role. Moreover, since 
quantum synchronization is imprinted in the spectral pu¬ 
rity of the emitted radiation [28] , the generated light may 
potentially serve as a direct diagnostic tool of quantum 
coherences in generic systems beyond cold gases such as 
organic molecules. 
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APPENDIX A: Mean-field approach 

The mean-field ansatz, p = 0 a p a , reduces the dy¬ 
namics to 3 N coupled nonlinear differential equations 
presented in the main text. In the most generic case 
we define local order parameters to take into account the 
effect of the inhomogeneous couplings: 

X a e^ = y]/(r a 6 )S'ft L e J06 , Y a e = '^2g(Y a b)Ste l(l>b . 

b^a b^a 

If the local order parameters vary slowly over the sys¬ 
tem size, and can be approximated to be the same 
for all dipoles one can define X a « / e ffZ, Y a « 
g e ffZ, where the global order parameter Z is defined as 
Ze 1 ® = S^e %< ^ a and the effective couplings are 

given by / e ff = Ea Eb^a f( r ab)/(N - 1 ) and g eff = 

E a E b*a9(Tab)/(N*»l). 

The steady-state solution Z = 0, <f> = Cot leads to two 
self-consistent equations for the order parameter Z and 
the collective frequency uJ 

ZP[f e ffQ + 2# e ff(J a + a;)] 

V N Ql< S * + ^) 2 + ( 2 /e ff ^ 2 + 2 ^eV 2 +Q 2 ]’ ' 

0= y^_ ZP[g e fiQ - 2 /eff(Jq + aJ)] _ 

V N QW« + ^) 2 + ( 2 /eff^ 2 + 2(&Z 2 + Q 2 Y } 

which can be evaluated in the N -> oo limit as integrals 
when the detunnings 5 a have a known distribution. 
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FIG. Al. Pair-wise two-time correlation functions in 
the steady-state parametrized by — Aexp(—qr). 

The correlations are calculated for a pair of dipoles in an en¬ 
semble of N — 70 dipoles identically coupled with /(r a &) fi§ 
fes/N. The cumulant expansion solution agrees with the full 
calculation except at W <C T where subradiant behavior dom¬ 
inates (purple region) [28]. 


APPENDIX B: Cumulant expansion approach and 
two-time correlation between dipoles 


The cumulant expansion method is a useful theoretical 
tool for including correlation effects beyond the mean- 
field approximation [28, 71, 72]. We keep two-point cor¬ 
relations such as i but f ac !° r i ze three-point 

correlations and higher [73]: 

(K°b°c) = (K^b)( a 2) + (^a)(^b a 2) + (< 3 'a<7)(< 5 f> 

-2 (a“)<af)(a?). (Bl) 


This factorization closes the set of dynamical equations 
of motion for all single particle observables and 

equal time two-point correlations. Two-time correlation 
functions can be computed by solving [48]: 


d(vj(t + T)a b (t)) 
dr 


iS a + 


r + B"i 


7 » (t + T)a b (t)) 


+ (* + T ) & b (*)) 

j^a,b 


where we have introduced the approximation (a*(t + 
r)a+ (t + r)a b (t)) « (^(f)}(<r+ (t + r)a b (t)). Compar- 
isons with exact numerical solutions show that the cumu¬ 
lant expansion captures well the steady-state behavior 
for inhomogeneous couplings /(r a &), provided the elas¬ 
tic couplings g(r a ij) are sufficiently small. In Fig. Al 
we compute the pair-wise two-time correlation function, 
z a,b{r) =lim t-,oo((o-+(t+r)+o-^(i+r))(o--(i)+CT^(f))), 
using both the cumulant expansion and the exact solu¬ 
tion. The decay rate of these correlations, Z a ^{r) = 
Ae -r7 , encodes information about the spectral coher¬ 
ence of the emitted radiation (note that here v = 0). 
The result shows that r/q exhibits the same dependence 
on W/Y as Zq. 


APPENDIX C: Conditional evolution, 
entanglement and quantum correlations 

An individual experimental realization can be consid¬ 
ered as a single trajectory, whose evolution can be quite 
different from the ensemble averaged solution of the mas¬ 
ter equation. Tracking the evolution of an individual 
trajectory is equivalent to performing continuous mea¬ 
surements that collect the record of the emitted photons, 
for example homodyne measurements. The conditional 
evolution of the system subject to continuous measure¬ 
ments can be modeled by the method of quantum state 
diffusion [48, 74]. For a single run the state of the system 
remains pure, p c = \ip) (^|, but the average over many 
trials reduces the system into a mixed state and recovers 
the density matrix obtained from the master equation. 

To probe the entanglement of the dipoles, in Fig. 2 
we calculate the average quantum Fisher information for 
each individual trajectory [46, 47] 

Fq(Pc) = ~^(Fq(p c ; J x ) + Fq(/ 5 c ; J y ) + Fg(p c ; J z )), 

where F Q (p c ;H) = 4((V>| U 2 |f/>) - (ip\ U |f/>) 2 ), and J XtVtZ 
are collective angular momentum operators. 

States with zero entanglement can still be nonclas- 
sical. Two systems are correlated if they share infor¬ 
mation with each other. The total amount of correla¬ 
tion can be quantified by the quantum mutual informa¬ 
tion X = Sa + Sb — Sab , where Si is the von Neu¬ 
mann entropy of the subsystem i G { A , T>, AB} (AB 
is the total system spanned by A and B together) , 
Si = — Tr[pilog 2 Pi], with pi the reduced density matrix 
of the subsystem i. A value varying between 0 and 2 
is obtained when A and B are pure states or maximally 
correlated respectively. The mutual information can be 
separated into a classical and a quantum part. The clas¬ 
sical part is Jb\a = max{<Se — Sb\a}- Here Sb\a is 
the von Neumann entropy of subsystem B conditioned 
on the measurement performed on A and max represents 
maximum value obtainable over all local measurements 
on A. The quantum part, known as the quantum dis¬ 
cord, T>b\a —F — Jb\a-> measures the amount of corre¬ 
lations that exceed the classical part and characterizes 
the “quantumness” of the system [55]. A state with 
nonzero quantum discord behaves in a way intrinsically 
non-classical, since a local measurement performed on 
one of its subsystems can disturb the whole system. In or¬ 
der to calculate Jb\a-> we consider a set of von Neumann 
measurements n ^ =1 2 = |(1 d= • & A ) with |r7 ^| 2 = 1 , 
made on the subsystem A and minimize the correspond- 
ing conditional entropy S b \a = Efc=i Tr bfc 5 '(/ 5 B|n£)], 
where p k = Tr[fl£p], p B \n£ = Tr A[n^p]/p fe [55]. In 
Fig. 2(e) we calculate the mutual information from X and 
the quantum discord from V using as subsystems A and 
B a pair of dipoles, a and b respectively. 
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